*! version 1.0, 11-11-2011
capture program drop nddfdual
program define nddfdual,rclass
    version 16

    gettoken word 0 : 0, parse(" =:,")
    while `"`word'"' != ":" & `"`word'"' != "=" {
        if `"`word'"' == "," | `"`word'"'=="" {
                error 198
        }
        local invars `invars' `word'
        gettoken word 0 : 0, parse("=:,")
    }
    unab invars : `invars'

    gettoken word 0 : 0, parse(" =:,")
    while `"`word'"' != ":" & `"`word'"' != "=" {
        if `"`word'"' == "," | `"`word'"'=="" {
                error 198
        }
        local gopvars `gopvars' `word'
        gettoken word 0 : 0, parse(" =:,")
    }
    unab gopvars : `gopvars'
	
    syntax varlist [if] [in], dmu(varname) [Time(varname) SEQuential  SAVing(string) maxiter(numlist integer >0 max=1) tol(numlist max=1) wmat(string) NOCheck] 
   * check whether the new version is available
   if "`nocheck'"==""{ //default
      if "$nddfdual"==""   deadualcheckupdate nddfdual 
   }       
		
	marksample touse 
		

	local bopvars `varlist'
	
	local invars: list uniq invars
	local gopvars: list uniq gopvars
	local bopvars: list uniq bopvars
	
	confirm numeric var `invars' `gopvars' `bopvars'
	
	local comvars: list invars & gopvars 
	if !(`"`comvars'"'==""){
		disp as error "`comvars' should not be specified as input and desriable output simultaneously."
		error 498
	}
	
	local comvars: list invars & bopvars
	if !(`"`comvars'"'==""){
		disp as error "`comvars' should not be specified as input and undesriable output simultaneously."
		error 498
	}	
	
	local comvars: list gopvars & bopvars
	if !(`"`comvars'"'==""){
		disp as error "`comvars' should not be specified as desriable and undesriable outputs simultaneously."
		error 498
	}	
		
  
	local techflag "=="
  
	if "`sequential'"!=""{
		if "`time'"==""{
		   disp as error "For sequential SBM model, time() should be specified."
		   error 498
		}
		else{
		   local techflag "<="
		}
	
	}
	
	
	
    preserve
	qui keep   `invars' `gopvars' `bopvars' `dmu' `time' `touse'
	qui gen _Row=_n
	label var _Row "Row #"
	qui keep if `touse'		
	
	foreach v in `invars' `gopvars' `bopvars'{
	   qui gen double p_`v'=.
	   label var p_`v' `"virtual price:`v'"'
	   local slackvars `slackvars' p_`v'
	}	
	
	
	
	tempvar rflag dmu2 tvar
	
	if `"`time'"'==""{
		
		qui gen `rflag'=1
		_nddfdual if `touse', rflag(`rflag') gen(`slackvars') invars(`invars') gopvars(`gopvars') bopvars(`bopvars') maxiter(`maxiter') tol(`tol') wmat(`wmat')
	}
	else{
	    qui egen `tvar'=group(`time')
		qui egen `dmu2'=group(`dmu')
		
		sort `dmu2' `tvar' _Row
		
		qui su `tvar',meanonly
		
		local tmax=r(max)
		qui gen `rflag'=0
		forv j=1/`tmax'{
			qui replace `rflag'=1 if `tvar' `techflag' `j'
			
			_nddfdual if `touse' & `tvar'==`j', rflag(`rflag') gen(`slackvars') invars(`invars') gopvars(`gopvars') bopvars(`bopvars') maxiter(`maxiter') tol(`tol')  wmat(`wmat')
			
		}
		
	}
	
    order _Row `dmu' `time'  `slackvars'
	keep _Row `dmu' `time'  `slackvars'
	
	disp _n(2) " Nddf dual LPs' results:"
	disp "    (_Row: Row # in the original data)"
	//disp "      S_X : Slack of X"
	list _Row `dmu' `time' `slackvars', sep(0) 

	//disp _n
	if `"`saving'"'!=""{
	  save `saving'
	  gettoken filenames saving:saving, parse(",")
	  local filenames `filenames'.dta
	  disp _n `"Estimated Results are saved in `filenames'."'
	}
	//tempname resmat
	//mkmat _Row `dmu' `time' TE `slackvars', mat(`resmat')
	//matrix list `resmat', noblank nohalf  noheader f(%9.6g)
	//return mat results=`resmat'
	return local file `filenames'
	restore 
	

end


	





capture program drop _nddfdual
program define _nddfdual
    version 16

    syntax  [if] [in], [rflag(varname) maxiter(numlist integer >0 max=1) tol(numlist max=1)  wmat(string)] gen(string) invars(varlist) gopvars(varlist) bopvars(varlist)
        marksample touse 
		markout `touse' `invars' `gopvars' `bopvars'
		
		tempvar touse2
		mark `touse2' if `rflag'
		markout `touse2' `invars' `gopvars' `bopvars'
		//qui gen `touse2'=`rflag'	

        local data `invars' `gopvars' `bopvars'
        local k: word count `invars'
        local l: word count `gopvars'
		local q: word count `bopvars'
		local nvar = `k'+`l'+`q'
		if "`maxiter'"==""{
			local maxiter=-1
		}
		if "`tol'"==""{
			local tol=-1
		}
		
	tempname weightvec
	
	if `"`wmat'"'!=""{
		confirm matrix `wmat'
		local ncol=colsof(`wmat')
		if `ncol'!=`nvar'{
		    dis as error `"# of column of `wmat' != # of input-output variables"'
			exit 498
		}
		mat `weightvec'=`wmat'
		
		
	   }
	else{
		
		mat  `weightvec'=(J(1,`k',1)*(1/3/`k'),J(1,`l',1)*(1/3/`l'),J(1,`q',1)*(1/3/`q'))
			
	}
		
		
		
		mata: wmat=st_matrix("`weightvec'")
        mata: nddfdual2("`data'","`touse'", "`touse2'",`k',`l',"`gen'",`maxiter',`tol',wmat')

end   



cap mata mata drop lpres()
mata:
    struct lpres   { real scalar    fval 
                     real matrix    coeff
                     real scalar    converged
                     real scalar    returncode
                      }
end


cap mata mata drop nddfdual()
mata:
   function nddfdual(real colvector    X, ///
                    real colvector    Y,  ///
                    real colvector    B, ///
                    real matrix    Xref, ///
                    real matrix    Yref, ///
                    real matrix    Bref, ///
					real scalar    maxiter, ///
					real scalar     tol , ///
					real colvector W)
   {
	
        class LinearProgram scalar q
		
		gx=-X
		gy=Y
		gb=-B
        nx=length(X)
        ny=length(Y)
        nb=length(B)
		wx=W[1::nx]
		wy=W[(nx+1)::(nx+ny)]
		wb=W[(nx+ny+1)::(nx+ny+nb)]
        s=ny+nb
        c=(X',-Y',B')
        lowerbd=J(1,length(c),0)
        upperbd=J(1,length(c),.)
	
        q = LinearProgram()
		q.setMaxOrMin("min")
        //A1 = (Xref',-Yref',Bref')
		A1 = (-Xref',Yref',-Bref')
        b1 = J(cols(Xref),1,0)
        A2 = (diag(gx),J(nx,ny,0),J(nx,nb,0))
		
        b2 =-wx
        A3 = (J(ny,nx,0),-diag(gy),J(ny,nb,0))
        b3 =-wy
		
        A4 = (J(nb,nx,0),J(nb,ny,0),diag(gb))
        b4 =-wb                   
        Aie=A1 \ A2 \A3 \A4
        bie=b1 \ b2 \b3 \ b4
        q.setCoefficients(c)
        q.setInequality(Aie, bie)
        q.setBounds(lowerbd, upperbd)
		
		if(maxiter!=-1){
		  q.setMaxiter(maxiter)
		}
		if (tol!=-1){
		  q.setTol(tol)
		}		
		
        theta=q.optimize()
        
         struct lpres scalar retres
         retres.fval=q.value()
         retres.coeff=q.parameters()
         retres.converged=q.converged()
         retres.returncode=q.returncode()
         return(retres)

}

end


cap mata mata drop nddfdual2()
mata:
void function nddfdual2(string scalar d, ///
                       string scalar touse, ///
                       string scalar rflag, ///
                       real scalar k,  ///
                       real scalar l,  ///
                       string scalar name, ///
					   real scalar maxiter, ///
					   real scalar  tol, ///
					   real colvector W)
    { 

          data=st_data(.,d,touse)
          data=data'
          dataref=st_data(.,d,rflag)
          dataref=dataref'
          M=rows(data)
          Xref=dataref[1..k,.]
          Yref=dataref[k+1..k+l,.]
          Bref=dataref[k+l+1..M,.]
          X=data[1..k,.]
          Y=data[k+1..k+l,.]
          B=data[k+l+1..M,.]		  

          beta=J(cols(data),M,.)
          struct lpres scalar res
          for(j=1;j<=cols(data);j++){
            res=nddfdual(X[.,j],Y[.,j],B[.,j],Xref,Yref,Bref,maxiter,tol,W)
            if(res.converged==1){
                beta[j,.]=res.coeff
              }
            }

          st_view(gen=.,.,name,touse)
          gen[.,.]=beta
    
    }

end

